From: tsteven4 <13596209+tsteven4@users.noreply.github.com> Date: Wed, 7 Aug 2024 21:04:13 +0000 (-0600) Subject: unify access to ellipsoid parameters (#1312) X-Git-Tag: archive/raspbian/1.10.0+ds-2+rpi1~1^2~12^2^2~78 X-Git-Url: https://dgit.raspbian.org/%22http:/www.example.com/cgi/%22https://%22Program/%22http:/www.example.com/cgi/%22https:/%22Program?a=commitdiff_plain;h=3230ea5e6c7abf48c37bedb3833ff95d74d9666c;p=gpsbabel.git unify access to ellipsoid parameters (#1312) * define constants for often used spherioid parms. * add function to calculate semi-minor axis. * use semi major axis function more. * have semi_minor_axis consume an ellipse * compute semi minor axis in GPS_Ellipse class. add constants for other ellipses that are used directly. --- diff --git a/jeeps/gpsdatum.h b/jeeps/gpsdatum.h index 2cfbb366a..eea7a3612 100644 --- a/jeeps/gpsdatum.h +++ b/jeeps/gpsdatum.h @@ -6,14 +6,33 @@ struct GPS_Ellipse { const char* name; double a; double invf; + + constexpr double b() const { + return a - a/invf; + } }; +// EPSG:7001 +constexpr GPS_Ellipse Airy1830_Ellipse = { "Airy 1830", 6377563.396, 299.3249646 }; + +// EPSG:7002 +constexpr GPS_Ellipse Airy1830Modified_Ellipse = { "Airy 1830 Modified", 6377340.189, 299.3249646 }; + +// EPSG:7004 +constexpr GPS_Ellipse Bessel1841_Ellipse = { "Bessel 1841", 6377397.155, 299.1528128 }; + +// EPSG:7019 +constexpr GPS_Ellipse GRS80_Ellipse = { "GRS80", 6378137.000, 298.257222101 }; + +// EPSG:4326 +constexpr GPS_Ellipse WGS84_Ellipse = { "WGS84", 6378137.000, 298.257223563 }; + const GPS_Ellipse GPS_Ellipses[]= { - { "Airy 1830", 6377563.396, 299.3249646 }, - { "Airy 1830 Modified", 6377340.189, 299.3249646 }, + Airy1830_Ellipse, + Airy1830Modified_Ellipse, { "Australian National", 6378160.000, 298.25 }, { "Bessel 1841 (Namibia)", 6377483.865, 299.1528128 }, - { "Bessel 1841", 6377397.155, 299.1528128 }, + Bessel1841_Ellipse, { "Clarke 1866", 6378206.400, 294.9786982 }, { "Clarke 1880", 6378249.145, 293.465 }, { "Everest (India 1830)", 6377276.345, 300.8017 }, @@ -30,12 +49,12 @@ const GPS_Ellipse GPS_Ellipses[]= { { "Krassovsky 1940", 6378245.000, 298.3 }, { "GRS67", 6378160.000, 6356774.516 }, { "GRS75", 6378140.000, 6356755.288 }, - { "GRS80", 6378137.000, 298.257222101 }, + GRS80_Ellipse, { "S. American 1969", 6378160.000, 298.25 }, { "WGS60", 6378165.000, 298.3 }, { "WGS66", 6378145.000, 298.25 }, { "WGS72", 6378135.000, 298.26 }, - { "WGS84", 6378137.000, 298.257223563 }, + WGS84_Ellipse, { "Clarke 1880 (Benoit)", 6378300.789, 293.466 }, }; diff --git a/jeeps/gpsmath.cc b/jeeps/gpsmath.cc index c97b6ad48..6edfb5a14 100644 --- a/jeeps/gpsmath.cc +++ b/jeeps/gpsmath.cc @@ -429,8 +429,8 @@ void GPS_Math_XYZ_To_LatLonH(double* phi, double* lambda, double* H, void GPS_Math_Airy1830LatLonH_To_XYZ(double phi, double lambda, double H, double* x, double* y, double* z) { - double a = 6377563.396; - double b = 6356256.910; + constexpr double a = Airy1830_Ellipse.a; + constexpr double b = Airy1830_Ellipse.b(); GPS_Math_LatLonH_To_XYZ(phi,lambda,H,x,y,z,a,b); @@ -456,8 +456,8 @@ void GPS_Math_Airy1830LatLonH_To_XYZ(double phi, double lambda, double H, void GPS_Math_WGS84LatLonH_To_XYZ(double phi, double lambda, double H, double* x, double* y, double* z) { - double a = 6378137.000; - double b = 6356752.314245; + constexpr double a = WGS84_Ellipse.a; + constexpr double b = WGS84_Ellipse.b(); GPS_Math_LatLonH_To_XYZ(phi,lambda,H,x,y,z,a,b); @@ -483,8 +483,8 @@ void GPS_Math_WGS84LatLonH_To_XYZ(double phi, double lambda, double H, void GPS_Math_XYZ_To_Airy1830LatLonH(double* phi, double* lambda, double* H, double x, double y, double z) { - double a = 6377563.396; - double b = 6356256.910; + constexpr double a = Airy1830_Ellipse.a; + constexpr double b = Airy1830_Ellipse.b(); GPS_Math_XYZ_To_LatLonH(phi,lambda,H,x,y,z,a,b); @@ -509,8 +509,8 @@ void GPS_Math_XYZ_To_Airy1830LatLonH(double* phi, double* lambda, double* H, void GPS_Math_XYZ_To_WGS84LatLonH(double* phi, double* lambda, double* H, double x, double y, double z) { - double a = 6378137.000; - double b = 6356752.314245; + constexpr double a = WGS84_Ellipse.a; + constexpr double b = WGS84_Ellipse.b(); GPS_Math_XYZ_To_LatLonH(phi,lambda,H,x,y,z,a,b); @@ -643,8 +643,8 @@ void GPS_Math_Airy1830M_LatLonToINGEN(double phi, double lambda, double* E, double F0 = 1.000035; double phi0 = 53.5; double lambda0 = -8.; - double a = 6377340.189; - double b = 6356034.447; + constexpr double a = Airy1830Modified_Ellipse.a; + constexpr double b = Airy1830Modified_Ellipse.b(); GPS_Math_LatLon_To_EN(E,N,phi,lambda,N0,E0,phi0,lambda0,F0,a,b); @@ -674,8 +674,8 @@ void GPS_Math_Airy1830LatLonToNGEN(double phi, double lambda, double* E, double F0 = 0.9996012717; double phi0 = 49.; double lambda0 = -2.; - double a = 6377563.396; - double b = 6356256.910; + constexpr double a = Airy1830_Ellipse.a; + constexpr double b = Airy1830_Ellipse.b(); GPS_Math_LatLon_To_EN(E,N,phi,lambda,N0,E0,phi0,lambda0,F0,a,b); @@ -704,7 +704,7 @@ int32_t GPS_Math_WGS84_To_Swiss_EN(double lat, double lon, double* E, const double lambda0 = 7.43958333; const double E0 = 600000.0; const double N0 = 200000.0; - double phi, lambda, alt, a, b; + double phi, lambda, alt; if (lat < 44.89022757) { return 0; @@ -713,9 +713,8 @@ int32_t GPS_Math_WGS84_To_Swiss_EN(double lat, double lon, double* E, return 0; } - assert(strcmp(GPS_Ellipses[4].name, "Bessel 1841") == 0); - a = GPS_Ellipses[4].a; - b = a - (a / GPS_Ellipses[4].invf); + constexpr double a = Bessel1841_Ellipse.a; + constexpr double b = Bessel1841_Ellipse.b(); GPS_Math_WGS84_To_Known_Datum_M(lat, lon, 0, &phi, &lambda, &alt, 123); GPS_Math_Swiss_LatLon_To_EN(phi, lambda, E, N, phi0, lambda0, E0, N0, a, b); @@ -742,11 +741,10 @@ void GPS_Math_Swiss_EN_To_WGS84(double E, double N, double* lat, double* lon) const double lambda0 = 7.43958333; const double E0 = 600000.0; const double N0 = 200000.0; - double phi, lambda, alt, a, b; + double phi, lambda, alt; - assert(strcmp(GPS_Ellipses[4].name, "Bessel 1841") == 0); - a = GPS_Ellipses[4].a; - b = a - (a / GPS_Ellipses[4].invf); + constexpr double a = Bessel1841_Ellipse.a; + constexpr double b = Bessel1841_Ellipse.b(); GPS_Math_Swiss_EN_To_LatLon(E, N, &phi, &lambda, phi0, lambda0, E0, N0, a, b); GPS_Math_Known_Datum_To_WGS84_M(phi, lambda, 0, lat, lon, &alt, 123); @@ -1112,7 +1110,7 @@ int32_t GPS_Math_WGS84_To_ICS_EN(double lat, double lon, double* E, int32_t ellipse = GPS_Datums[datum].ellipse; a = GPS_Ellipses[ellipse].a; - b = a - (a / GPS_Ellipses[ellipse].invf); + b = GPS_Ellipses[ellipse].b(); GPS_Math_WGS84_To_Known_Datum_M(lat, lon, 0, &phi, &lambda, &alt, datum); GPS_Math_Cassini_LatLon_To_EN(phi, lambda, E, N, @@ -1148,7 +1146,7 @@ void GPS_Math_ICS_EN_To_WGS84(double E, double N, double* lat, double* lon) int32_t ellipse = GPS_Datums[datum].ellipse; a = GPS_Ellipses[ellipse].a; - b = a - (a / GPS_Ellipses[ellipse].invf); + b = GPS_Ellipses[ellipse].b(); GPS_Math_Cassini_EN_To_LatLon(E, N, &phi, &lambda, phi0, lambda0, E0, N0, a, b); @@ -1307,8 +1305,8 @@ void GPS_Math_NGENToAiry1830LatLon(double E, double N, double* phi, double F0 = 0.9996012717; double phi0 = 49.; double lambda0 = -2.; - double a = 6377563.396; - double b = 6356256.910; + constexpr double a = Airy1830_Ellipse.a; + constexpr double b = Airy1830_Ellipse.b(); GPS_Math_EN_To_LatLon(E,N,phi,lambda,N0,E0,phi0,lambda0,F0,a,b); @@ -1337,8 +1335,8 @@ void GPS_Math_INGENToAiry1830MLatLon(double E, double N, double* phi, double F0 = 1.000035; double phi0 = 53.5; double lambda0 = -8.; - double a = 6377340.189; - double b = 6356034.447; + constexpr double a = Airy1830Modified_Ellipse.a; + constexpr double b = Airy1830Modified_Ellipse.b(); GPS_Math_EN_To_LatLon(E,N,phi,lambda,N0,E0,phi0,lambda0,F0,a,b); @@ -1539,15 +1537,13 @@ void GPS_Math_Known_Datum_To_WGS84_M(double Sphi, double Slam, double SH, { double Sa; double Sif; - double Da; - double Dif; double x; double y; double z; int32_t idx; - Da = 6378137.0; - Dif = 298.257223563; + constexpr double Da = WGS84_Ellipse.a; + constexpr double Dif = WGS84_Ellipse.invf; idx = GPS_Datums[n].ellipse; Sa = GPS_Ellipses[idx].a; @@ -1581,8 +1577,6 @@ void GPS_Math_WGS84_To_Known_Datum_M(double Sphi, double Slam, double SH, double* Dphi, double* Dlam, double* DH, int32_t n) { - double Sa; - double Sif; double Da; double Dif; double x; @@ -1590,8 +1584,8 @@ void GPS_Math_WGS84_To_Known_Datum_M(double Sphi, double Slam, double SH, double z; int32_t idx; - Sa = 6378137.0; - Sif = 298.257223563; + constexpr double Sa = WGS84_Ellipse.a; + constexpr double Sif = WGS84_Ellipse.invf; idx = GPS_Datums[n].ellipse; Da = GPS_Ellipses[idx].a; @@ -1628,9 +1622,6 @@ void GPS_Math_Known_Datum_To_WGS84_C(double Sphi, double Slam, double SH, double Sa; double Sif; double Sb; - double Da; - double Dif; - double Db; double x; double y; double z; @@ -1639,9 +1630,8 @@ void GPS_Math_Known_Datum_To_WGS84_C(double Sphi, double Slam, double SH, double sy; double sz; - Da = 6378137.0; - Dif = 298.257223563; - Db = Da - (Da / Dif); + constexpr double Da = WGS84_Ellipse.a; + constexpr double Db = WGS84_Ellipse.b(); idx = GPS_Datums[n].ellipse; Sa = GPS_Ellipses[idx].a; @@ -1682,23 +1672,19 @@ void GPS_Math_WGS84_To_Known_Datum_C(double Sphi, double Slam, double SH, double* Dphi, double* Dlam, double* DH, int32_t n) { - double Sa; - double Sif; double Da; double Dif; double x; double y; double z; int32_t idx; - double Sb; double Db; double dx; double dy; double dz; - Sa = 6378137.0; - Sif = 298.257223563; - Sb = Sa - (Sa / Sif); + constexpr double Sa = WGS84_Ellipse.a; + constexpr double Sb = WGS84_Ellipse.b(); idx = GPS_Datums[n].ellipse; Da = GPS_Ellipses[idx].a; @@ -1803,9 +1789,7 @@ void GPS_Math_Known_Datum_To_Known_Datum_C(double Sphi, double Slam, double SH, double* DH, int32_t n1, int32_t n2) { double Sa; - double Sif; double Da; - double Dif; double x1; double y1; double z1; @@ -1824,8 +1808,7 @@ void GPS_Math_Known_Datum_To_Known_Datum_C(double Sphi, double Slam, double SH, idx1 = GPS_Datums[n1].ellipse; Sa = GPS_Ellipses[idx1].a; - Sif = GPS_Ellipses[idx1].invf; - Sb = Sa - (Sa / Sif); + Sb = GPS_Ellipses[idx1].b(); x1 = GPS_Datums[n1].dx; y1 = GPS_Datums[n1].dy; @@ -1833,8 +1816,7 @@ void GPS_Math_Known_Datum_To_Known_Datum_C(double Sphi, double Slam, double SH, idx2 = GPS_Datums[n2].ellipse; Da = GPS_Ellipses[idx2].a; - Dif = GPS_Ellipses[idx2].invf; - Db = Da - (Da / Dif); + Db = GPS_Ellipses[idx2].b(); x2 = GPS_Datums[n2].dx; y2 = GPS_Datums[n2].dy; @@ -2116,8 +2098,6 @@ int32_t GPS_Math_NAD83_To_UTM_EN(double lat, double lon, double* E, double N0; double E0; double F0; - double a; - double b; if (!GPS_Math_LatLon_To_UTM_Param(lat,lon,zone,zc,&lambda0,&E0, &N0,&F0)) { @@ -2126,9 +2106,8 @@ int32_t GPS_Math_NAD83_To_UTM_EN(double lat, double lon, double* E, phi0 = 0.0; - assert(strcmp(GPS_Ellipses[21].name, "GRS80") == 0); - a = GPS_Ellipses[21].a; - b = a - (a / GPS_Ellipses[21].invf); + constexpr double a = GRS80_Ellipse.a; + constexpr double b = GRS80_Ellipse.b(); GPS_Math_LatLon_To_EN(E,N,lat,lon,N0,E0,phi0,lambda0,F0,a,b); @@ -2294,7 +2273,7 @@ int32_t GPS_Math_Known_Datum_To_UTM_EN(double lat, double lon, double* E, idx = GPS_Datums[n].ellipse; a = GPS_Ellipses[idx].a; - b = a - (a / GPS_Ellipses[idx].invf); + b = GPS_Ellipses[idx].b(); GPS_Math_LatLon_To_EN(E,N,lat,lon,N0,E0,phi0,lambda0,F0,a,b); @@ -2517,7 +2496,7 @@ void GPS_Math_UTM_EN_to_LatLon(int ReferenceEllipsoid, //found at http://www.gpsy.com/gpsinfo/geotoutm/index.html double k0 = 0.9996; - double a, b; + double a, f; double eccSquared; double eccPrimeSquared; double e1; @@ -2526,8 +2505,8 @@ void GPS_Math_UTM_EN_to_LatLon(int ReferenceEllipsoid, double x, y; a = GPS_Ellipses[ReferenceEllipsoid].a; - b = 1 / GPS_Ellipses[ReferenceEllipsoid].invf; - eccSquared = b * (2.0 - b); + f = 1 / GPS_Ellipses[ReferenceEllipsoid].invf; + eccSquared = f * (2.0 - f); e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared)); x = UTMEasting - E0; //remove false easting diff --git a/reference/grid-bng~csv.gpx b/reference/grid-bng~csv.gpx index 0c361caec..ddcb600ea 100644 --- a/reference/grid-bng~csv.gpx +++ b/reference/grid-bng~csv.gpx @@ -1,8 +1,8 @@ - - + + AlineLodge Aline Lodge @@ -16,7 +16,7 @@ Caernarfon City (Small) - + Camborne Camborne @@ -37,28 +37,28 @@ Forfar City (Small) - + Hawick Hawick Hawick City (Small) - + Hosta Hosta Hosta City (Small) - + IsleofRhum Isle of Rhum Isle of Rhum City (Small) - + Lerwick Lerwick